 % Simulation
 
clc
clear all
close all

% calibrated beta and d_I
beta=0.15151; % 0.1682 for low (80) / 0.15151 for high (95)
d_I=0.00768; % 0.03576 for low / 0.00768 for high

w=ones(26,2); % no wage differentials
load moments_from_Seoul_survey.mat
load distance_GIS.mat
load calibrated_E.mat
load initial_infection.mat

epsilon_weekday=4.1642;
epsilon_weekend=4.9144;
kappa=0.0339;

H_R(:,1)=H_R_young; % young
H_R(:,2)=H_R_old;   % old

ncity=25;

gamma=1/18; % 18 days to be recovered
tau=[1/8.5 1/10.2]; % days to be quaurantined
death_rate=[0.0021 0.0273]; % fatality rate

delta_t(:,1)=[0.00583956102 0.00583956102]'; % young, weekday weekend
delta_t(:,2)=[0.00736463186 0.00736463186]'; % old, weekday weekend

delta_j(:,1)=[0.00209103789 0.00138626867 0.00209103789 0.00138626867]'; % young, weekday case/visit weekend case/visit
delta_j(:,2)=[0.00247622592 0.00096650016 0.00247622592 0.00096650016]'; % old, weekday case/visit weekend case/visit
delta_j=zeros(4,2); % no disclosure

global w d_GIS ncity epsilon_weekday epsilon_weekend kappa E_weekday E_weekend H_R initial_infection gamma tau death_rate beta d_I delta_j delta_t

time=1000; period=[1:1:time];

%% Actual Simulation
% setup dates: starting from Jan 30
Day = datetime(2020,1,30) + caldays(0:time);

% Initial Values
    Qtime(1,:,1)=initial_infection'; % Initial detected, as of January 30 (all looks exogenous)
    % city 8(1), 13(1), 15(1), 16(1) // 4 the detected
    Qtime(1,:,2)=zeros(1,25);
    Itime(1,:,1)=9*initial_infection';
    Itime(1,:,2)=zeros(1,25);
    for a=1:2
        for d1=1:ncity
            Ntime(1,d1,a)=H_R(d1,a);
            Stime(1,d1,a)=Ntime(1,d1,a)-Itime(1,d1,a)-Qtime(1,d1,a);
            Rtime(1,d1,a)=Ntime(1,d1,a)-Itime(1,d1,a)-Stime(1,d1,a)-Qtime(1,d1,a); % zero
            Ftime(1,d1,a)=death_rate(a)*Rtime(1,d1,a);
            Dtime(1,d1,a)=0;
            cum_Dtime(1,d1)=Qtime(1,d1,1)+Qtime(1,d1,2);
            cum_Dtime_14(1,d1)=Qtime(1,d1,1)+Qtime(1,d1,2);
        end
    end

% run simulation
for t=2:time
    [DayNumber,DayName] = weekday(Day(t));
    if DayNumber>=2 && DayNumber<=6 % weekdays
            pi(t,:,:,:)=Calculate_pi_function(cum_Dtime_14(t-1,:));
                for d1=1:ncity
                    for a=1:2
                    S_dot(t,d1,a)=0;
                    for d2=1:ncity
                        S_dot(t,d1,a)=S_dot(t,d1,a)-beta*(sum(pi(t,:,d2,1).*Itime(t-1,:,1))+sum(pi(t,:,d2,2).*Itime(t-1,:,2)))/(sum(pi(t,:,d2,1).*(Ntime(t-1,:,1)-Qtime(t-1,:,1)))+sum(pi(t,:,d2,2).*(Ntime(t-1,:,2)-Qtime(t-1,:,2))))*(pi(t,d1,d2,a)*Stime(t-1,d1,a));
                    end
                    Stime(t,d1,a) =  Stime(t-1,d1,a) + S_dot(t,d1,a);
                    Itime(t,d1,a) =  Itime(t-1,d1,a)*(1 - gamma) - S_dot(t,d1,a) - d_I*Itime(t-1,d1,a);
                    Qtime(t,d1,a) =  Qtime(t-1,d1,a)*(1 - tau(a)) + d_I*Itime(t-1,d1,a);
                    Rtime(t,d1,a) =  Rtime(t-1,d1,a) + Itime(t-1,d1,a)*gamma + tau(a)*Qtime(t-1,d1,a);
                    Ftime(t,d1,a) =  death_rate(a)*Rtime(t,d1,a);
                    Ntime(t,d1,a) =  H_R(d1,a) - death_rate(a)*Rtime(t,d1,a);
                    Dtime(t,d1,a) = d_I*Itime(t-1,d1,a);
                    end
                    cum_Dtime(t,d1) = cum_Dtime(t-1,d1) + d_I*Itime(t-1,d1,1) + d_I*Itime(t-1,d1,2); % adding newly detected cases
                    cum_Dtime_14(t,d1) = sum(Dtime(max(1,t-13):t,d1,1))+sum(Dtime(max(1,t-13):t,d1,2));
                end
    else % weekends
            pi(t,:,:,:)=Calculate_pi_weekend_function(cum_Dtime_14(t-1,:));
                for d1=1:ncity
                    for a=1:2
                    S_dot(t,d1,a)=0;
                    for d2=1:ncity
                        S_dot(t,d1,a)=S_dot(t,d1,a)-beta*(sum(pi(t,:,d2,1).*Itime(t-1,:,1))+sum(pi(t,:,d2,2).*Itime(t-1,:,2)))/(sum(pi(t,:,d2,1).*(Ntime(t-1,:,1)-Qtime(t-1,:,1)))+sum(pi(t,:,d2,2).*(Ntime(t-1,:,2)-Qtime(t-1,:,2))))*(pi(t,d1,d2,a)*Stime(t-1,d1,a));
                    end
                    Stime(t,d1,a) =  Stime(t-1,d1,a) + S_dot(t,d1,a);
                    Itime(t,d1,a) =  Itime(t-1,d1,a)*(1 - gamma) - S_dot(t,d1,a) - d_I*Itime(t-1,d1,a);
                    Qtime(t,d1,a) =  Qtime(t-1,d1,a)*(1 - tau(a)) + d_I*Itime(t-1,d1,a);
                    Rtime(t,d1,a) =  Rtime(t-1,d1,a) + Itime(t-1,d1,a)*gamma + tau(a)*Qtime(t-1,d1,a);
                    Ftime(t,d1,a) =  death_rate(a)*Rtime(t,d1,a);
                    Ntime(t,d1,a) =  H_R(d1,a) - death_rate(a)*Rtime(t,d1,a);
                    Dtime(t,d1,a) = d_I*Itime(t-1,d1,a);
                    end
                    cum_Dtime(t,d1) = cum_Dtime(t-1,d1) + d_I*Itime(t-1,d1,1) + d_I*Itime(t-1,d1,2); % adding newly detected cases
                    cum_Dtime_14(t,d1) = sum(Dtime(max(1,t-13):t,d1,1))+sum(Dtime(max(1,t-13):t,d1,2));
                end
    end
end

cum_Dtime_in_two_month=cum_Dtime(123,:)';
target1=sum(cum_Dtime_in_two_month) % 861

Itime_sum=Itime(:,:,1)+Itime(:,:,2);
Qtime_sum=Qtime(:,:,1)+Qtime(:,:,2);
fraction=sum(Itime_sum,2)./(sum(Qtime_sum,2)+sum(Itime_sum,2)); % Our primary estimates for the fraction of infections that are undetected range from 88.7% to 93.6%. (Iceland)
%[temp peak]=max(sum(Qtime_sum,2));
target2=mean(fraction(1:123)) % new target 90%

day=period;
Pop=sum(H_R);
Stime_sum=Stime(:,:,1)+Stime(:,:,2);
Rtime_sum=Rtime(:,:,1)+Rtime(:,:,2);
Ftime_sum=Ftime(:,:,1)+Ftime(:,:,2);

common_z_weekday=draw_frechet(E_weekday,epsilon_weekday,10^5);
common_z_weekend=draw_frechet(E_weekend,epsilon_weekend,10^5);

global common_z_weekday common_z_weekend
for t=1:time
    [DayNumber,DayName] = weekday(Day(t));
    if DayNumber>=2 && DayNumber<=6 % weekdays
        utility_weekdays(t,:,:)=Calculate_utility_per_capita_function(cum_Dtime_14(t,:),squeeze(Qtime(t,:,:)),squeeze(Rtime(t,:,:)));
    else
        utility_weekends(t,:,:)=Calculate_utility_per_capita_weekend_function(cum_Dtime_14(t,:),squeeze(Qtime(t,:,:)),squeeze(Rtime(t,:,:)));
    end
    t
end

utility_weekdays_sum(:,1)=sum(utility_weekdays(:,:,1),2); % per capita
utility_weekdays_sum(:,2)=sum(utility_weekdays(:,:,2),2); % per capita
size_weekdays=size(utility_weekdays_sum);
for a=1:2
for t=1:size_weekdays(1)
    if utility_weekdays_sum(t,a)==0
        utility_weekdays_sum(t,a)=NaN;
    end
end
end
utility_weekdays_sum_norm(:,1)=log(utility_weekdays_sum(:,1)./utility_weekdays_sum(1,1));
utility_weekdays_sum_norm(:,2)=log(utility_weekdays_sum(:,2)./utility_weekdays_sum(1,2));

utility_weekends_sum(:,1)=sum(utility_weekends(:,:,1),2); % per capita
utility_weekends_sum(:,2)=sum(utility_weekends(:,:,2),2); % per capita
size_weekends=size(utility_weekends_sum);
for a=1:2
for t=1:size_weekends(1)
    if utility_weekends_sum(t,a)==0
        utility_weekends_sum(t,a)=NaN;
    end
end
end
utility_weekends_sum_norm(:,1)=log(utility_weekends_sum(:,1)./utility_weekends_sum(3,1)); % first weekend, day 3
utility_weekends_sum_norm(:,2)=log(utility_weekends_sum(:,2)./utility_weekends_sum(3,2)); % first weekend, day 3

utility_weekdays_sum_norm_avg=utility_weekdays_sum_norm(:,1)*(Pop(1)/(Pop(1)+Pop(2)))+utility_weekdays_sum_norm(:,2)*(Pop(2)/(Pop(1)+Pop(2)));
utility_weekends_sum_norm_avg=utility_weekends_sum_norm(:,1)*(Pop(1)/(Pop(1)+Pop(2)))+utility_weekends_sum_norm(:,2)*(Pop(2)/(Pop(1)+Pop(2)));

for t=2:time
    [DayNumber,DayName] = weekday(Day(t));
    if DayNumber>=2 && DayNumber<=6 % weekdays
        utility_sum_norm_avg(t,1)=utility_weekdays_sum_norm_avg(t);
        utility_sum_norm_young(t,1)=utility_weekdays_sum_norm(t,1);
        utility_sum_norm_old(t,1)=utility_weekdays_sum_norm(t,2);
    else
        utility_sum_norm_avg(t,1)=utility_weekends_sum_norm_avg(t);
        utility_sum_norm_young(t,1)=utility_weekends_sum_norm(t,1);
        utility_sum_norm_old(t,1)=utility_weekends_sum_norm(t,2);
    end
end

%% Save files
total_cum_Dtime=sum(cum_Dtime,2);
%save ('Output/sen_detection_low_moment_original.mat', 'day', 'Pop', 'Stime', 'Itime', 'Qtime', 'Rtime', 'Ftime', 'total_cum_Dtime', 'utility_sum_norm_avg', 'utility_sum_norm_young','utility_sum_norm_old')
%save ('Output/sen_detection_low_moment_delta1_0_delta2_0.mat', 'day', 'Pop', 'Stime', 'Itime', 'Qtime', 'Rtime', 'Ftime', 'total_cum_Dtime', 'utility_sum_norm_avg', 'utility_sum_norm_young','utility_sum_norm_old')
%save ('Output/sen_detection_high_moment_original.mat', 'day', 'Pop', 'Stime', 'Itime', 'Qtime', 'Rtime', 'Ftime', 'total_cum_Dtime', 'utility_sum_norm_avg', 'utility_sum_norm_young','utility_sum_norm_old')
save ('Output/sen_detection_high_moment_delta1_0_delta2_0.mat', 'day', 'Pop', 'Stime', 'Itime', 'Qtime', 'Rtime', 'Ftime', 'total_cum_Dtime', 'utility_sum_norm_avg', 'utility_sum_norm_young','utility_sum_norm_old')